! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! +                                                           +
! +  glint_timestep.f90 - part of the Glimmer-CISM ice model  + 
! +                                                           +
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! 
! Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009, 2010
! Glimmer-CISM contributors - see AUTHORS file for list of contributors
!
! This file is part of Glimmer-CISM.
!
! Glimmer-CISM is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 2 of the License, or (at
! your option) any later version.
!
! Glimmer-CISM is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with Glimmer-CISM.  If not, see <http://www.gnu.org/licenses/>.
!
! Glimmer-CISM is hosted on BerliOS.de:
! https://developer.berlios.de/projects/glimmer-cism/
!
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#ifdef HAVE_CONFIG_H
#include "config.inc"
#endif

module glint_timestep
  !*FD timestep of a GLINT instance

  use glint_type
  use glint_constants
  private
  public glint_i_tstep

contains

  subroutine glint_i_tstep(time,            instance,       &
                           g_temp,          g_temp_range,   &
                           g_precip,        g_zonwind,      &
                           g_merwind,       g_humid,        &
                           g_lwdown,        g_swdown,       &
                           g_airpress,      g_orog,         &
                           g_orog_out,      g_albedo,       &
                           g_ice_frac,      g_veg_frac,     &
                           g_snowice_frac,  g_snowveg_frac, &
                           g_snow_depth,                    &
                           g_water_in,      g_water_out,    &
                           t_win,           t_wout,         &
                           ice_vol,         out_f,          &
                           orogflag,        ice_tstep,      &
                           gcm_smb_in,                      &
                           qsmb_g,          tsfc_g,         &
                           topo_g,          gmask,          &
                           gfrac,           gtopo,          &
                           grofi,           grofl,          &
                           ghflx )

    !*FD Performs time-step of an ice model instance. 
    !*FD Note that input quantities here are accumulated/average totals since the
    !*FD last call.
    !
    ! TO DO - The interface requires many arguments that may not be needed.
    ! Some of these arguments (e.g., zonwind, merwind, humid, lwdown, swdown, airpress)
    ! should be made optional.
    
    use glide
    use glide_setup
    use glide_io
    use glimmer_paramets
    use glint_io
    use glint_mbal_io
    use glint_climate
    use glimmer_routing
    use glimmer_log
    use glimmer_physcon, only: rhow,rhoi
    use glide_mask, only: glide_mask_ocean

    implicit none

    ! ------------------------------------------------------------------------  
    ! Arguments
    ! ------------------------------------------------------------------------  

    integer,                intent(in)   :: time         !*FD Current time in hours
    type(glint_instance), intent(inout)  :: instance     !*FD Model instance
    real(rk),dimension(:,:),intent(in)   :: g_temp       !*FD Global mean surface temperature field ($^{\circ}$C)
    real(rk),dimension(:,:),intent(in)   :: g_temp_range !*FD Global surface temperature half-range field ($^{\circ}$C)
    real(rk),dimension(:,:),intent(in)   :: g_precip     !*FD Global precip field total (mm)

    real(rk),dimension(:,:),intent(in)   :: g_zonwind    !*FD Global mean surface zonal wind (m/s)
    real(rk),dimension(:,:),intent(in)   :: g_merwind    !*FD Global mean surface meridonal wind (m/s)
    real(rk),dimension(:,:),intent(in)   :: g_humid      !*FD Global surface humidity (%)
    real(rk),dimension(:,:),intent(in)   :: g_lwdown     !*FD Global downwelling longwave (W/m^2)
    real(rk),dimension(:,:),intent(in)   :: g_swdown     !*FD Global downwelling shortwave (W/m^2)
    real(rk),dimension(:,:),intent(in)   :: g_airpress   !*FD Global surface air pressure (Pa)

    real(rk),dimension(:,:),intent(in)   :: g_orog       !*FD Input global orography (m)
    real(rk),dimension(:,:),intent(out)  :: g_orog_out   !*FD Output orography (m)
    real(rk),dimension(:,:),intent(out)  :: g_albedo     !*FD Output surface albedo 
    real(rk),dimension(:,:),intent(out)  :: g_ice_frac   !*FD Output ice fraction
    real(rk),dimension(:,:),intent(out)  :: g_veg_frac   !*FD Output veg fraction
    real(rk),dimension(:,:),intent(out)  :: g_snowice_frac !*FD Output snow-ice fraction
    real(rk),dimension(:,:),intent(out)  :: g_snowveg_frac !*FD Output snow-veg fraction
    real(rk),dimension(:,:),intent(out)  :: g_snow_depth !*FD Output snow depth (m)
    real(rk),dimension(:,:),intent(out)  :: g_water_in   !*FD Input water flux (m)
    real(rk),dimension(:,:),intent(out)  :: g_water_out  !*FD Output water flux (m)
    real(rk),               intent(out)  :: t_win        !*FD Total water input (kg)
    real(rk),               intent(out)  :: t_wout       !*FD Total water output (kg)
    real(rk),               intent(out)  :: ice_vol      !*FD Output ice volume (m$^3$)
    type(output_flags),     intent(in)   :: out_f        !*FD Flags to tell us whether to do output   
    logical,                intent(in)   :: orogflag     !*FD Set if we have new global orog
    logical,                intent(out)  :: ice_tstep    !*FD Set if we have done an ice time step
    logical,                  optional,intent(in)  :: gcm_smb_in ! true if getting sfc mass balance from a GCM
    real(rk),dimension(:,:,:),optional,intent(in)  :: qsmb_g    ! Depth of new ice (m)
    real(rk),dimension(:,:,:),optional,intent(in)  :: tsfc_g    ! Surface temperature (C)
    real(rk),dimension(:,:,:),optional,intent(in)  :: topo_g    ! Surface elevation (m)
    integer, dimension(:,:),  optional,intent(in)  :: gmask     ! = 1 where global data are valid, else = 0
    real(rk),dimension(:,:,:),optional,intent(out) :: gfrac     ! ice fractional area [0,1]
    real(rk),dimension(:,:,:),optional,intent(out) :: gtopo     ! surface elevation (m)
    real(rk),dimension(:,:,:),optional,intent(out) :: grofi     ! ice runoff (kg/m^2/s = mm H2O/s)
    real(rk),dimension(:,:,:),optional,intent(out) :: grofl     ! liquid runoff (kg/m^2/s = mm H2O/s)
    real(rk),dimension(:,:,:),optional,intent(out) :: ghflx     ! heat flux (W/m^2, positive down)

    ! ------------------------------------------------------------------------  
    ! Internal variables
    ! ------------------------------------------------------------------------  

    real(rk),dimension(:,:),pointer :: upscale_temp => null() ! temporary array for upscaling
    real(rk),dimension(:,:),pointer :: routing_temp => null() ! temporary array for flow routing
    real(rk),dimension(:,:),pointer :: accum_temp   => null() ! temporary array for accumulation
    real(rk),dimension(:,:),pointer :: ablat_temp   => null() ! temporary array for ablation
    integer, dimension(:,:),pointer :: fudge_mask   => null() ! temporary array for fudging
    real(sp),dimension(:,:),pointer :: thck_temp    => null() ! temporary array for volume calcs
    real(sp),dimension(:,:),pointer :: calve_temp   => null() ! temporary array for calving flux
    real(rk) :: start_volume,end_volume,flux_fudge
    integer :: i
    integer :: j, ii, jj, nx, ny, il, jl

    logical :: gcm_smb   ! true if getting sfc mass balance from a GCM

    gcm_smb = .false.
    if (present(gcm_smb_in)) then
       gcm_smb = gcm_smb_in
    endif

    ! Zero outputs

    if (present(gfrac)) gfrac(:,:,:) = 0._rk
    if (present(gtopo)) gtopo(:,:,:) = 0._rk
    if (present(grofi)) grofi(:,:,:) = 0._rk
    if (present(grofl)) grofl(:,:,:) = 0._rk
    if (present(ghflx)) ghflx(:,:,:) = 0._rk

    if (GLC_DEBUG) then
       write(stdout,*) ' '
       write(stdout,*) 'In glint_i_tstep, time =', time
       write(stdout,*) 'next_time =', instance%next_time
    endif

    ! Check whether we're doing anything this time.

    if (time /= instance%next_time) then
       return
    else
       instance%next_time = instance%next_time + instance%mbal_tstep
    end if

    ! Assume we always need this, as it's too complicated to work out when we do and don't

    call coordsystem_allocate(instance%lgrid, thck_temp)
    call coordsystem_allocate(instance%lgrid, calve_temp)
    ice_tstep = .false.

    ! Downscale input fields from global to local grid
    ! This subroutine computes instance%acab and instance%artm, the key inputs to GLIDE.

    if (gcm_smb) then
       call glint_downscaling_gcm (instance,              &
                                   qsmb_g,      tsfc_g,   &
                                   topo_g,      gmask)
    else
       call glint_downscaling(instance,                  &
                              g_temp,     g_temp_range,  &
                              g_precip,   g_orog,        &
                              g_zonwind,  g_merwind,     &
                              g_humid,    g_lwdown,      &
                              g_swdown,   g_airpress,    &
                              orogflag)
    endif

    ! ------------------------------------------------------------------------  
    ! Sort out some local orography and remove bathymetry. This relies on the 
    ! point 1,1 being underwater. However, it's a better method than just 
    ! setting all points < 0.0 to zero
    ! ------------------------------------------------------------------------  

    call glide_get_usurf(instance%model, instance%local_orog)
    call glint_remove_bath(instance%local_orog,1,1)


    if (.not. gcm_smb) then

       ! ------------------------------------------------------------------------  
       ! Adjust the surface temperatures using the lapse-rate, by reducing to
       ! sea-level and then back up to high-res orography
       ! ------------------------------------------------------------------------  

       call glint_lapserate(instance%artm, real(instance%global_orog,rk), real(-instance%data_lapse_rate,rk))
       call glint_lapserate(instance%artm, real(instance%local_orog,rk),  real(instance%lapse_rate,rk))

       ! Process the precipitation field if necessary ---------------------------
       ! and convert from mm/s to m/s

       call glint_calc_precip(instance)

    endif   ! not gcm_smb

    ! Get ice thickness ----------------------------------------

    call glide_get_thk(instance%model,thck_temp)

    ! Do accumulation --------------------------------------------------------

    if (gcm_smb) then
       call glint_accumulate_gcm(instance%mbal_accum, time, instance%acab, instance%artm)
    else
       call glint_accumulate(instance%mbal_accum, time, instance%artm, instance%arng, instance%prcp, &
                             instance%snowd, instance%siced, instance%xwind, instance%ywind, &
                             instance%local_orog, real(thck_temp,rk), instance%humid,    &
                             instance%swdown, instance%lwdown, instance%airpress)
    endif

    ! Initialise water budget quantities to zero. These will be over-ridden if
    ! there's an ice-model time-step

    t_win=0.0       ; t_wout=0.0
    g_water_out=0.0 ; g_water_in=0.0

    if (GLC_DEBUG) then
       write(stdout,*) ' '
       write(stdout,*) 'Check for ice dynamics timestep'
       write(stdout,*) 'time =', time
       write(stdout,*) 'start_time =', instance%mbal_accum%start_time
       write(stdout,*) 'mbal_step =', instance%mbal_tstep
       write(stdout,*) 'mbal_accum_time =', instance%mbal_accum_time
    endif

    ! ------------------------------------------------------------------------  
    ! ICE TIMESTEP begins HERE ***********************************************
    ! ------------------------------------------------------------------------  

    if (time - instance%mbal_accum%start_time + instance%mbal_tstep == instance%mbal_accum_time) then

       if (instance%mbal_accum_time < instance%ice_tstep) then 
          instance%next_time = instance%next_time + instance%ice_tstep - instance%mbal_tstep
       end if

       ice_tstep = .true.

       ! Prepare arrays for water budgeting

       if (out_f%water_out .or. out_f%total_wout .or. out_f%water_in .or. out_f%total_win) then
          call coordsystem_allocate(instance%lgrid, accum_temp)
          call coordsystem_allocate(instance%lgrid, ablat_temp)
          accum_temp = 0.0
          ablat_temp = 0.0
       end if

       ! Calculate the initial ice volume (scaled and converted to water equivalent)
       call glide_get_thk(instance%model, thck_temp)
       thck_temp = thck_temp*real(rhoi/rhow)
       start_volume = sum(thck_temp)

       ! ---------------------------------------------------------------------
       ! do the different parts of the glint timestep
       ! ---------------------------------------------------------------------

       do i = 1, instance%n_icetstep

          if (GLC_DEBUG) then
             write (stdout,*) 'GLIDE timestep, iteration =', i
          endif

          ! Calculate the initial ice volume (scaled and converted to water equivalent)
          call glide_get_thk(instance%model,thck_temp)
          thck_temp = thck_temp*real(rhoi/rhow)

          ! Get latest upper-surface elevation (needed for masking)
          call glide_get_usurf(instance%model, instance%local_orog)
          call glint_remove_bath(instance%local_orog,1,1)

          ! Get the mass-balance, as m water/year 
          call glint_get_mbal(instance%mbal_accum, instance%artm, instance%prcp, instance%ablt, &
                              instance%acab, instance%snowd, instance%siced, instance%mbal_accum_time)

          ! Mask out non-accumulation in ice-free areas

          where(thck_temp<=0.0 .and. instance%acab<0.0)
             instance%acab = 0.0
             instance%ablt = instance%prcp
          end where

          ! Set acab to zero for ocean cells (bed below sea level, no ice present)

          where (instance%model%geometry%thkmask == glide_mask_ocean)
             instance%acab = 0.0
          endwhere

          ! Put climate inputs in the appropriate places, with conversion ----------

          ! Note on units: 
          ! For this subroutine, input acab is in m/yr; this value is multiplied 
          !  by tim0/(scyr*thk0) and copied to data%climate%acab.
          ! Input artm is in deg C; this value is copied to data%climate%artm (no unit conversion).

          call glide_set_acab(instance%model, instance%acab*real(rhow/rhoi))
          call glide_set_artm(instance%model, instance%artm)

          if (GLC_DEBUG) then
             il = itest_local
             jl = jtest_local
             write (stdout,*) ' '
             write (stdout,*) 'After glide_set_acab, glide_set_artm: i, j =', il, jl
             write (stdout,*) 'acab (m/y), artm (C) =', instance%acab(il,jl), instance%artm(il,jl)
          endif

          ! Adjust glint acab and ablt for output
 
          where (instance%acab < -thck_temp .and. thck_temp > 0.0)
             instance%acab = -thck_temp
             instance%ablt =  thck_temp
          end where

          instance%glide_time=instance%glide_time+instance%model%numerics%tinc
          call glide_tstep_p1(instance%model,instance%glide_time)
          call glide_tstep_p2(instance%model)
          call glide_tstep_p3(instance%model,no_write=.true.)

          ! Add the calved ice to the ablation field

!lipscomb - TO DO - Calving needs to be added to the ice runoff.
!                   Also add basal melting (bmlt) to the liquid runoff.

          call glide_get_calving(instance%model, calve_temp)
          calve_temp = calve_temp * real(rhoi/rhow)

          instance%ablt = instance%ablt + calve_temp/instance%model%numerics%tinc
          instance%acab = instance%acab - calve_temp/instance%model%numerics%tinc

          ! Accumulate for water-budgeting
          if (out_f%water_out .or. out_f%total_wout .or. out_f%water_in .or. out_f%total_win) then
             accum_temp = accum_temp + instance%prcp*instance%model%numerics%tinc
             ablat_temp = ablat_temp + instance%ablt*instance%model%numerics%tinc
          endif

          ! Save GLIDE output until now
          call glide_io_writeall(instance%model,instance%model)
          call glint_io_writeall(instance,instance%model)

       end do

       ! Calculate flux fudge factor --------------------------------------------

       if (out_f%water_out .or. out_f%total_wout .or. out_f%water_in .or. out_f%total_win) then

          call coordsystem_allocate(instance%lgrid,fudge_mask)

          call glide_get_thk(instance%model,thck_temp)
          end_volume=sum(thck_temp)

          where (thck_temp>0.0)
             fudge_mask = 1
          elsewhere
             fudge_mask = 0
          endwhere

          flux_fudge = (start_volume + sum(accum_temp) - sum(ablat_temp) - end_volume) / sum(fudge_mask)

          ! Apply fudge_factor

          where(thck_temp > 0.0)
             ablat_temp = ablat_temp + flux_fudge
          endwhere
          
          deallocate(fudge_mask)
          fudge_mask => null()

       endif

       ! Upscale water flux fields ----------------------------------------------
       ! First water input (i.e. mass balance + ablation)

       if (out_f%water_in) then
          call coordsystem_allocate(instance%lgrid, upscale_temp)

          where (thck_temp > 0.0)
             upscale_temp = accum_temp
          elsewhere
             upscale_temp = 0.0
          endwhere

          call mean_to_global(instance%ups,   &
                              upscale_temp,   &
                              g_water_in,     &
                              instance%out_mask)
          deallocate(upscale_temp)
          upscale_temp => null()
       endif

       ! Now water output (i.e. ablation) - and do routing

!lipscomb - TO DO - Modify (or skip?) the following for gcm_smb option

       if (out_f%water_out) then
          call coordsystem_allocate(instance%lgrid, upscale_temp)
          call coordsystem_allocate(instance%lgrid, routing_temp)

          where (thck_temp > 0.0)
             upscale_temp = ablat_temp
          elsewhere
             upscale_temp = 0.0
          endwhere

          call glide_get_usurf(instance%model, instance%local_orog)
          call flow_router(instance%local_orog, &
                           upscale_temp, &
                           routing_temp, &
                           instance%out_mask, &
                           real(instance%lgrid%delta%pt(1),rk), &
                           real(instance%lgrid%delta%pt(2),rk))

          call mean_to_global(instance%ups,   &
                              routing_temp,   &
                              g_water_out,    &
                              instance%out_mask)

          deallocate(upscale_temp,routing_temp)
          upscale_temp => null()
          routing_temp => null()

       endif

       ! Sum water fluxes and convert if necessary ------------------------------

       if (out_f%total_win) then
          t_win  = sum(accum_temp) * instance%lgrid%delta%pt(1)* &
                                     instance%lgrid%delta%pt(2)
       endif

       if (out_f%total_wout) then
          t_wout = sum(ablat_temp) * instance%lgrid%delta%pt(1)* &
                                     instance%lgrid%delta%pt(2)
       endif

    end if

    ! Output instantaneous values

    call glint_mbal_io_writeall(instance, instance%model,       &
                                outfiles = instance%out_first,  &
                                time = real(time*hours2years,sp))

    ! ------------------------------------------------------------------------ 
    ! Upscaling of output
    ! ------------------------------------------------------------------------ 

    ! We now upscale all fields at once...

    call get_i_upscaled_fields(instance, g_orog_out, g_albedo, g_ice_frac, g_veg_frac, &
                               g_snowice_frac, g_snowveg_frac, g_snow_depth)

    ! Calculate ice volume ---------------------------------------------------

    if (out_f%ice_vol) then
       call glide_get_thk(instance%model, thck_temp)
       ice_vol = sum(thck_temp) * instance%lgrid%delta%pt(1)* &
                                  instance%lgrid%delta%pt(2)
    endif

    ! Tidy up ----------------------------------------------------------------

    if (associated(accum_temp)) then 
       deallocate(accum_temp)
       accum_temp => null()
    end if

    if (associated(ablat_temp)) then
       deallocate(ablat_temp)
       ablat_temp => null()
    end if

    if (associated(calve_temp)) then
       deallocate(calve_temp)
       calve_temp => null()
    end if

    deallocate(thck_temp)
    thck_temp => null()

  end subroutine glint_i_tstep

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine glint_remove_bath(orog,x,y)

    !*FD Sets ocean areas to zero height, working recursively from
    !*FD a known ocean point.

    real(sp),dimension(:,:),intent(inout) :: orog !*FD Orography --- used for input and output
    integer,                intent(in)    :: x,y  !*FD Location of starting point (index)

    integer :: nx,ny

    nx=size(orog,1) ; ny=size(orog,2)

    if (orog(x,y).lt.0.0) orog(x,y)=0.0
    call glint_find_bath(orog,x,y,nx,ny)

  end subroutine glint_remove_bath

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  recursive subroutine glint_find_bath(orog,x,y,nx,ny)

    !*FD Recursive subroutine called by {\tt glimmer\_remove\_bath}.

    real(sp),dimension(:,:),intent(inout) :: orog  !*FD Orography --- used for input and output
    integer,                intent(in)    :: x,y   !*FD Starting point
    integer,                intent(in)    :: nx,ny !*FD Size of array {\tt orography}

    integer,dimension(4) :: xi=(/ -1,1,0,0 /)
    integer,dimension(4) :: yi=(/ 0,0,-1,1 /)
    integer :: ns=4,i

    do i=1,ns
       if (x+xi(i).le.nx.and.x+xi(i).gt.0.and. &
            y+yi(i).le.ny.and.y+yi(i).gt.0) then
          if (orog(x+xi(i),y+yi(i)).lt.0.0) then
             orog(x+xi(i),y+yi(i))=0.0
             call glint_find_bath(orog,x+xi(i),y+yi(i),nx,ny)
          endif
       endif
    enddo

  end subroutine glint_find_bath

end module glint_timestep


